Detecting Seasonality in population indicators
#load(paste0(IO$output_kindara,"pop_indicators.Rdata"), verbose = TRUE)
load(paste0(IO$tmp_kindara, "sk1.Rdata"), verbose = TRUE)
## Loading objects:
## sk1
g = ggplot(sk1, aes(x = date, y = tot_sex)) + geom_line()
g
s = reshape(sk1, varying = list(4:ncol(sk1)), idvar = c("date","n_cycles","n_obs"), direction = "long")
s$time = colnames(sk1)[4:ncol(sk1)][s$time]
colnames(s) = c("date","n_cycles","n_obs","indicator","value")
g = ggplot(s[grep("sex",s$indicator),], aes(x = date, y = value/n_cycles, col = indicator))
g + geom_line() + ggtitle("Sexual frequency (#/number of cycles)")
g = ggplot(s[grep("sex",s$indicator),], aes(x = date, y = value/n_obs, col = indicator))
g + geom_line() + ggtitle("Sexual frequency (#/number of observations on that day)")
s$n_indicator = s$value/s$n_obs
s$n_indicator = s$value/s$n_cycles
signal = s$n_indicator[s$indicator == "tot_sex"]
time = par$time_num - 1
plot(time, signal, type = "l")
seasonal_analysis = detect_seasonal_pattern(time = time, signal = signal, title = "tot_sex Kindara", remove_weekly_patterns = FALSE)
seasonal_analysis = detect_seasonal_pattern(time = time, signal = signal, title = "tot_sex Kindara", remove_weekly_patterns = TRUE)
seasonal_analysis$time_series$date = par$date_seq - 365
save(seasonal_analysis, file = paste0(IO$tmp_kindara, "seasonal_analysis_tot_sex.Rdata"))
### TRYING WITH STL (Seasonal Decomposition of Time Series)
## weekly + yearly
# weekly period
signal_ts = ts(signal, frequency = 7)
stl_w = stl(signal_ts, s.window = "periodic")
plot(stl_w)
# yearly period
signal_ts_no_weekly = ts(apply((stl_w$time.series[,2:3]),1,sum), frequency = length(signal)/3)
stl_y = stl(signal_ts_no_weekly, s.window = "periodic")
plot(stl_y)
# bi-annual period
signal_ts_no_weekly = ts(apply((stl_w$time.series[,2:3]),1,sum), frequency = length(signal)/6)
stl_ba = stl(signal_ts_no_weekly, s.window = "periodic")
plot(stl_ba)
## "detrend" + weekly + yearly
# yearly period >> to detrend
signal_ts = ts(signal, frequency = length(signal)/3)
stl_d = stl(signal_ts, s.window = "periodic")
plot(stl_d)
# weekly period
signal_detrended = ts(apply((stl_d$time.series[,2:3]),1,sum), frequency = 7)
stl_w = stl(signal_detrended, s.window = "periodic")
plot(stl_w)
# yearly period
signal_ts_no_weekly = ts(apply((stl_w$time.series[,2:3]),1,sum), frequency = length(signal)/3)
stl_y = stl(signal_ts_no_weekly, s.window = "periodic")
plot(stl_y)
### TRYING WITH Decompose
signal_ts = ts(signal, frequency = length(signal)/3)
decompose_y = decompose(signal_ts)
plot(decompose_y)
# yearly
signal_ts = ts(signal, frequency = length(signal)/3)
plot(signal_ts)
y = findfrequency(signal_ts)
y
## [1] 1
# weekly
signal_ts = ts(signal, frequency = 7)
plot(signal_ts)
y = findfrequency(signal_ts)
y
## [1] 1
# gaussian noise (negative example check)
signal_ts = ts(rnorm(length(signal)), frequency = 7)
plot(signal_ts)
y = findfrequency(signal_ts)
y
## [1] 1
# cosine of period 7 (positive example check)
signal_ts = ts(cos((1:length(signal))*2*pi/7), frequency = 7)
plot(signal_ts)
y = findfrequency(signal_ts)
y
## [1] 7
# cosine of period 7 (positive example check with noise)
signal_ts = ts(cos((1:length(signal))*2*pi/7)+rnorm(length(signal))/2, frequency = 7)
plot(signal_ts)
y = findfrequency(signal_ts)
y
## [1] 7
# trying to remove the trend first to check that it's not just the first component that predominates
seasonal_analysis = detect_seasonal_pattern(time = time, signal = signal, title = "tot_sex Kindara", remove_weekly_patterns = TRUE)
detrended_signal = signal - seasonal_analysis$time_series$trend
# yearly
detrended_signal_ts = ts(detrended_signal, frequency = length(signal)/3)
plot(detrended_signal_ts)
y = findfrequency(detrended_signal_ts)
y
## [1] 1
# weekly
detrended_signal_ts = ts(detrended_signal, frequency = 7)
plot(detrended_signal_ts)
y = findfrequency(detrended_signal_ts)
y
## [1] 1
# does not work either
fft_sex = fft(signal)
plot(abs(fft_sex), type = "l")
fft_sex = fft(detrended_signal)
plot(abs(fft_sex), type = "l")
period_in_days = length(signal)/(1:(length(signal)/2))
plot(period_in_days,abs(fft_sex)[2:(length(signal)/2+1)], type = "l", log = "x")
